
rm(list = ls())
### libraries and packages
pkg <- c("plyr", "dplyr", "tidyr", 
         "MASS", "multiwayvcov", "fracdiff",
         "fractal","compactr",
         "lme4", "ArfimaMLM","gam",
         "forecast",
         "stargazer", "lmtest", "doBy", "ggplot2", 
         "mediation", "vars", "DataCombine", "foreign",
         "mgcv", "ordinal", "reshape2", 
         "xtable", "sandwich", "rms")
lapply(pkg, require, character.only = TRUE)


# rid of E
options(scipen=999)
### setting wd for inputs
setwd("~/Dropbox/Pollution and Public Perceptions in China/Data and Analysis/Data")
load("ppp_cleaned")

# marginal effect for l_watch: 

model1 <- lm(l_watch ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov1 <- cluster.vcov(model1, p2$day) 
coeftest(model1, m.vcov1) 

beta.hat1 <- coef(model1) 

xi1 <- model.matrix(model1)[,2] # aqi used in the model

x01 <- seq(min(xi1), max(xi1), length.out = 1000) 

dy.dx1 <- beta.hat1["aqi"] + 2*beta.hat1["I(aqi^2)"]*x01

se.dy.dx1 <- sqrt(m.vcov1["aqi", "aqi"] + 4*(x01^2)*m.vcov1["I(aqi^2)", "I(aqi^2)"] + 4*x01*m.vcov1["aqi", "I(aqi^2)"])

upr1 <- dy.dx1 + 1.64*se.dy.dx1 # upper bound of the 90 percent 
lwr1 <- dy.dx1 - 1.64*se.dy.dx1 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x01), ylim = mm(c(upr1, lwr1)), main = "Effect of pollution on demand for local oversight",
      ylab = "Marginal Effect of AQI", 
      xlab = "Standardized AQI")
lines(x01, dy.dx1, lwd = 3)
lines(x01, lwr1, lty = 3)
lines(x01, upr1, lty = 3)


model2 <- lm(c_watch ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov2 <- cluster.vcov(model2, p2$day)
coeftest(model2, m.vcov2) 
beta.hat2 <- coef(model2) 

xi2 <- model.matrix(model2)[,2] 


x02 <- seq(min(xi2), max(xi2), length.out = 1000) 

dy.dx2 <- beta.hat2["aqi"] + 2*beta.hat2["I(aqi^2)"]*x02

se.dy.dx2 <- sqrt(m.vcov2["aqi", "aqi"] + 4*(x02^2)*m.vcov2["I(aqi^2)", "I(aqi^2)"] + 4*x02*m.vcov2["aqi", "I(aqi^2)"])

upr2 <- dy.dx2 + 1.64*se.dy.dx2 # upper bound of the 90 percent 
lwr2 <- dy.dx2 - 1.64*se.dy.dx2 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x02), ylim = mm(c(upr2, lwr2)), 
      ylab = "Marginal Effect of AQI", main = "Effect of pollution on demand for central oversight",
      xlab = "Standardized AQI")
lines(x02, dy.dx2, lwd = 3)
lines(x02, lwr2, lty = 3)
lines(x02, upr2, lty = 3)


model3 <- lm(l_sat ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov3 <- cluster.vcov(model3, p2$day) 
coeftest(model3, m.vcov3) 

beta.hat3 <- coef(model3) 

xi3 <- model.matrix(model3)[,2] 

x03 <- seq(min(xi3), max(xi3), length.out = 1000) 

dy.dx3 <- beta.hat3["aqi"] + 2*beta.hat3["I(aqi^2)"]*x03

se.dy.dx3 <- sqrt(m.vcov3["aqi", "aqi"] + 4*(x03^2)*m.vcov3["I(aqi^2)", "I(aqi^2)"] + 4*x03*m.vcov3["aqi", "I(aqi^2)"])

upr3 <- dy.dx3 + 1.64*se.dy.dx3 # upper bound of the 90 percent 
lwr3 <- dy.dx3 - 1.64*se.dy.dx3 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x03), ylim = mm(c(upr3, lwr3)), main = "Effect of pollution on satisfaction with local government",
      ylab = "Marginal Effect of AQI", 
      xlab = "Standardized AQI")
lines(x03, dy.dx3, lwd = 3)
lines(x03, lwr3, lty = 3)
lines(x03, upr3, lty = 3)

# marginal effect for c_sat: 

model4 <- lm(c_sat ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov4 <- cluster.vcov(model4, p2$day) 
coeftest(model4, m.vcov4) 

beta.hat4 <- coef(model4) 

xi4 <- model.matrix(model4)[,2] 

x04 <- seq(min(xi4), max(xi4), length.out = 1000) 

dy.dx4 <- beta.hat4["aqi"] + 2*beta.hat4["I(aqi^2)"]*x04

se.dy.dx4 <- sqrt(m.vcov4["aqi", "aqi"] + 4*(x04^2)*m.vcov4["I(aqi^2)", "I(aqi^2)"] + 4*x04*m.vcov4["aqi", "I(aqi^2)"])

upr4 <- dy.dx4 + 1.64*se.dy.dx4 # upper bound of the 90 percent 
lwr4 <- dy.dx4 - 1.64*se.dy.dx4 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x04), ylim = mm(c(upr4, lwr4)), 
      ylab = "Marginal Effect of AQI", main = "Effect of pollution on satisfaction with central government",
      xlab = "Standardized AQI")
lines(x04, dy.dx4, lwd = 3)
lines(x04, lwr4, lty = 3)
lines(x04, upr4, lty = 3)

dy.dx5 <- beta.hat5["aqi"] + 2*beta.hat5["I(aqi^2)"]*x05

se.dy.dx5 <- sqrt(m.vcov5["aqi", "aqi"] + 4*(x05^2)*m.vcov5["I(aqi^2)", "I(aqi^2)"] + 4*x05*m.vcov5["aqi", "I(aqi^2)"])

upr5 <- dy.dx5 + 1.64*se.dy.dx5 # upper bound of the 90 percent 
lwr5 <- dy.dx5 - 1.64*se.dy.dx5 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x05), ylim = mm(c(upr5, lwr5)), 
      ylab = "Marginal Effect of AQI", main = "Effect of pollution on confidence in the economy",
      xlab = "Standardized AQI")
lines(x05, dy.dx5, lwd = 3)
lines(x05, lwr5, lty = 3)
lines(x05, upr5, lty = 3)


model6 <- lm(anti_west ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov6 <- cluster.vcov(model6, p2$day) 
coeftest(model6, m.vcov6) 

beta.hat6 <- coef(model6)

xi6 <- model.matrix(model6)[,2] 

x06 <- seq(min(xi6), max(xi6), length.out = 1000) 

dy.dx6 <- beta.hat6["aqi"] + 2*beta.hat6["I(aqi^2)"]*x06

se.dy.dx6 <- sqrt(m.vcov6["aqi", "aqi"] + 4*(x06^2)*m.vcov6["I(aqi^2)", "I(aqi^2)"] + 4*x06*m.vcov6["aqi", "I(aqi^2)"])

upr6 <- dy.dx6 + 1.64*se.dy.dx6 # upper bound of the 90 percent 
lwr6 <- dy.dx6 - 1.64*se.dy.dx6 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x06), ylim = mm(c(upr6, lwr6)), 
      ylab = "Marginal Effect of AQI", main = "Effect of pollution on belief in globalization as a US conspiracy", 
      xlab = "Standardized AQI")
lines(x06, dy.dx6, lwd = 3)
lines(x06, lwr6, lty = 3)
lines(x06, upr6, lty = 3)

 
model7 <- lm(life_sat ~ aqi + I(aqi^2) + aqi.lagged + aqi.lagged2 + hukou 
            + insider + soe + educ 
            + kids + married + income +
              gender + ccp + parade + national + age, data=p2)
m.vcov7 <- cluster.vcov(model7, p2$day) 
coeftest(model7, m.vcov7) 

beta.hat7 <- coef(model7) 

xi7 <- model.matrix(model7)[,2] 

x07 <- seq(min(xi7), max(xi7), length.out = 1000) 


dy.dx7 <- beta.hat7["aqi"] + 2*beta.hat7["I(aqi^2)"]*x07

se.dy.dx7 <- sqrt(m.vcov7["aqi", "aqi"] + 4*(x07^2)*m.vcov7["I(aqi^2)", "I(aqi^2)"] + 4*x07*m.vcov7["aqi", "I(aqi^2)"])

upr7 <- dy.dx7 + 1.64*se.dy.dx7 # upper bound of the 90 percent 
lwr7 <- dy.dx7 - 1.64*se.dy.dx7 # lower bound of the 90 percent

# plot the marginal effect using package compactr
eplot(xlim = mm(x07), ylim = mm(c(upr7, lwr7)), 
      ylab = "Marginal Effect of AQI", main = "Effect of pollution on life satisfaction",
      xlab = "Standardized AQI")
lines(x07, dy.dx7, lwd = 3)
lines(x07, lwr7, lty = 3)
lines(x07, upr7, lty = 3)


## FIGURE A9
setwd("~/Dropbox/Pollution and Public Perceptions in China/Manuscript/Figures")
pdf("curv_lin1.pdf", height=4, width=4)
eplot(xlim = mm(x01), ylim = mm(c(upr1, lwr1)), main = "Effect on demand for local oversight",
      ylab = "Marginal Effect of AQI", 
      xlab = "Standardized AQI")
lines(x01, dy.dx1, lwd = 3)
lines(x01, lwr1, lty = 3)
lines(x01, upr1, lty = 3)
dev.off()

pdf("curv_lin1_5.pdf", height=4, width=4)
eplot(xlim = mm(x02), ylim = mm(c(upr2, lwr2)), 
      ylab = "Marginal Effect of AQI", main = "Effect on demand for central oversight",
      xlab = "Standardized AQI")
lines(x02, dy.dx2, lwd = 3)
lines(x02, lwr2, lty = 3)
lines(x02, upr2, lty = 3)
dev.off()

pdf("curv_lin2.pdf", height=4, width=4)
eplot(xlim = mm(x03), ylim = mm(c(upr3, lwr3)), main = "Effect on satisfaction with local govt",
      ylab = "Marginal Effect of AQI", 
      xlab = "Standardized AQI")
lines(x03, dy.dx3, lwd = 3)
lines(x03, lwr3, lty = 3)
lines(x03, upr3, lty = 3)
dev.off()

pdf("curv_lin2_5.pdf", height = 4, width = 4)
eplot(xlim = mm(x04), ylim = mm(c(upr4, lwr4)), 
      ylab = "Marginal Effect of AQI", main = "Effect on satisfaction with central govt",
      xlab = "Standardized AQI")
lines(x04, dy.dx4, lwd = 3)
lines(x04, lwr4, lty = 3)
lines(x04, upr4, lty = 3)
dev.off()

pdf("curv_lin3_5.pdf", height = 4, width = 4)
eplot(xlim = mm(x06), ylim = mm(c(upr6, lwr6)), 
      ylab = "Marginal Effect of AQI", main = "Effect on belief in anti-west conspiracy", 
      xlab = "Standardized AQI")
lines(x06, dy.dx6, lwd = 3)
lines(x06, lwr6, lty = 3)
lines(x06, upr6, lty = 3)
dev.off()

pdf("curv_lin4.pdf", height = 4, width = 4)
eplot(xlim = mm(x07), ylim = mm(c(upr7, lwr7)), 
      ylab = "Marginal Effect of AQI", main = "Effect on life satisfaction",
      xlab = "Standardized AQI")
lines(x07, dy.dx7, lwd = 3)
lines(x07, lwr7, lty = 3)
lines(x07, upr7, lty = 3)
dev.off()

